## Chapter 3 (Shared Identities) 
## Our Common Bonds 
## Winter 2021 - 2022  

## Read in source file 
source("ocb_replication_file.r")

######################
## Open-Ended Codes ## 
######################

## Read in this data 
codes <- read_csv(file="data/open_ended_codes.csv") 
mean(codes$nonsense, na.rm=T) ## 7% say gibberish 

pdf(file="figures/chi_levendusky_fig03001.pdf",
    height = 11,
    width = 8.5)
codes %>% 
  filter(nonsense == 0) %>% 
  summarise(Freedom = 100*mean(freedom,na.rm=T),
            Opportunity = 100*mean(opportunity, na.rm=T),
            Government = 100*mean(govt, na.rm=T),
            Diversity = 100*mean(diversity, na.rm=T),
            Economy = 100*mean(economy, na.rm=T),
            Nature = 100*mean(nature, na.rm=T),
            Patriotism = 100*mean(patriotism, na.rm=T),
            Unclear = 100*mean(unclear, na.rm=T),
            `Not Proud` = 100*mean(reject,na.rm=T)) %>% 
  pivot_longer(everything(),
               names_to = "type",
               values_to = "percentage") %>% 
  ggplot(aes(x = reorder(type, -percentage), 
             y = percentage)) + 
  geom_bar(stat = "identity") + 
  xlab("Rationale") + 
  ylab("Percentage") + 
  theme_bw() + 
  theme(panel.grid = element_blank()) + 
  ggtitle("") 
dev.off()

#######################################
## 2015 American Identity Experiment ##
#######################################

## Read in the data from the 2015 American Identity Experiment   
american <- read_dta(file="data/jop_exp1_data_for_book.dta") %>% 
  clean_names() 

## Out-Party Feeling Thermometer 
pdf(file="figures/chi_levendusky_fig03002.pdf")
ggplot(data = american) + 
  geom_density(aes(x=outparty_ft, 
                   y=..count..,
                   group=as.factor(treatment), 
                   linetype=as.factor(treatment),
                   fill=as.factor(treatment)), alpha=0.35) + 
  theme_bw() + 
  ylab("") +
  xlab("Out-Party Feeling Thermometer") + 
  ggtitle("") + 
  scale_fill_grey(name="", 
                  labels=c("Control","Treatment"), 
                  start = 0.10, 
                  end = 0.8) + 
  scale_linetype(name="",
                 labels=c("Control","Treatment"),
                 limits=c(0,1)) + 
  theme(plot.title = element_text(hjust=0.5),
        legend.position = "bottom",
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid = element_blank()) 
dev.off() 

## show what effect that has on 0 degree ratings & ratings of 50 or higher 
american <- american %>%
  mutate(ft_zero = ifelse(outparty_ft == 0,1,0),
         ft_fifty = ifelse(outparty_ft > 49,1,0),
         less_ten = ifelse(outparty_ft < 11,1,0))

t.test(american$ft_zero ~ american$treatment) 
t.test(american$ft_fifty ~ american$treatment) 
t.test(american$less_ten ~ american$treatment)

## FT ratings of 0 go from 13% in the control to 9% in the treatment (25% relative decline) 
## FT ratings of 50+ go from 19% in the control to 28% in the treatment (nearly 50% increase!) 

## Intelligent 
pdf(file="figures/chi_levendusky_fig03003.pdf")
ggplot(data=american) + 
  geom_bar(aes(x=tr_intelligent, 
               y=100*..prop..,
               fill=as.factor(treatment)),
           position="dodge") + 
  xlab("") + 
  ylab("Percentage")  + 
  scale_x_continuous(breaks=c(1:5),
                     labels=c("Not at All \n Well","Not Too \n Well",
                              "Somewhat \n Well","Very \n Well","Extremely \n Well")) + 
  theme_bw() +   
  ggtitle("") + 
  scale_fill_grey(name="", 
                  labels=c("Control","Treatment"), 
                  start = 0.10, 
                  end = 0.8) + 
  theme(plot.title = element_text(hjust=0.5),
        legend.position = "bottom",
        panel.grid = element_blank()) 
dev.off() 

## American 
pdf(file="figures/chi_levendusky_fig03004.pdf")
ggplot(data=american) + 
  geom_bar(aes(x=tr_american, 
               y=100*..prop..,
               fill=as.factor(treatment)),
           position="dodge") + 
  xlab("") + 
  ylab("Percentage")  + 
  scale_x_continuous(breaks=c(1:5),
                     labels=c("Not at All \n Well","Not Too \n Well",
                              "Somewhat \n Well","Very \n Well","Extremely \n Well")) + 
  theme_bw() +   
  ggtitle("") + 
  scale_fill_grey(name="", 
                  labels=c("Control","Treatment"), 
                  start = 0.10, 
                  end = 0.8) + 
  theme(plot.title = element_text(hjust=0.5),
        legend.position = "bottom",
        panel.grid = element_blank()) 
dev.off() 

## Look at the relative values 
prop.table(table(american$tr_american[american$treatment == 0]))
prop.table(table(american$tr_american[american$treatment == 1]))
## Control: 8%, 18%, 30% 25%, 20% 
## Treatment: 5%, 13%, 31% 27%, 25% 
## big differences are the fall in not at all/not too well, adn the gain in extremely well 

## Figures 3.5 & 3.6 come from the NAES Data 
## see the file "chapter_three_naes.r" for those graphs 

#################################
## 2020 American Identity Data ##
#################################

## Read in the survey data  
pilot <- read_csv(file="data/american_identity_under_trump_public.csv") 
## read in the open-ended codes 
codes <- read_csv(file="data/coded_responses.csv") 

with_codes <- left_join(x = codes,
                        y = pilot,
                        by = "respondent_id") 
## should fully join 
anti_join(x = codes,
          y = pilot,
          by = "respondent_id") 

## Graph to mimic one for 2015 codes 
pdf(file = "figures/chi_levendusky_fig03007.pdf",
    height = 11,
    width = 8.5) 
codes %>% 
  filter(nonsense == 0) %>% 
  mutate(reject = neg_racism + neg_govt + neg_general) %>% 
  summarise(Freedom = 100*mean(freedom,na.rm=T),
            Opportunity = 100*mean(opportunity, na.rm=T),
            Government = 100*mean(government, na.rm=T),
            Diversity = 100*mean(diversity, na.rm=T),
            Economy = 100*mean(economy, na.rm=T),
            Nature = 100*mean(nature, na.rm=T),
            Patriotism = 100*mean(patriotism, na.rm=T),
            Unclear = 100*mean(unclear, na.rm=T),
            `Not Proud` = 100*mean(reject,na.rm=T)) %>% 
  pivot_longer(everything(),
               names_to = "type",
               values_to = "percentage") %>% 
  ggplot(aes(x = reorder(type, -percentage), 
             y = percentage)) + 
  geom_bar(stat = "identity") + 
  xlab("Rationale") + 
  ylab("Percentage") + 
  theme_bw() + 
  theme(panel.grid = element_blank()) + 
  ggtitle("") 
dev.off() 
## not being proud is now the #3 reason (~12%)! 

## Break down the "not proud" reasons 
codes %>% 
  filter(neg_racism == 1 | neg_govt == 1 | neg_general == 1) %>% 
  summarise(racism = 100*mean(neg_racism, na.rm=T),
            trump = 100*mean(neg_govt, na.rm=T),
            other = 100*mean(neg_general, na.rm=T))
## of those who reject, split pretty evenly between racism, Trump, and other factors 

## partisan differences in negative traits? 
with_codes <- with_codes %>% 
  mutate(reject_prompt = ifelse((neg_racism == 1 | neg_govt == 1 | neg_general == 1),1,0)) 
table(with_codes$party_id[with_codes$reject_prompt==1])
t.test(with_codes$white ~ with_codes$reject_prompt)
## 85% of those who reject the prompt are Democrats!  
## Whites are also less likely to reject the prompt 

## Code up the treatment 
pilot <- pilot %>% 
  mutate(control = ifelse(!is.na(q113),1,0), 
         american = ifelse(!is.na(q38),1,0),
         sept11 = ifelse(!is.na(q42),1,0), 
         treatment = ifelse(control==1,0,
                            ifelse(american==1,1,
                                   ifelse(sept11==1,2,NA))),
         white = ifelse(q19==1,1,0), 
         rep = ifelse(party_id >4,1,0), 
         white_rep = ifelse(white == 1 & rep == 1,1,0))

## code up the DVs 
pilot <- pilot %>% 
  mutate(other_therm = ifelse(!is.na(q45_2),q45_2,NA),  
         other_american = (-1*q50_1) + 6,
         other_intelligent = (-1*q50_2) + 6,
         other_honest = (-1*q50_3) + 6,
         other_openminded = (-1*q50_4) + 6,
         other_generous = (-1*q50_5) + 6,
         rev_bipart = (-1*q41) + 6) 

## do the traits scale? 
psych::alpha(cbind(pilot$other_american,
                   pilot$other_intelligent,
                   pilot$other_honest,
                   pilot$other_openminded,
                   pilot$other_generous,
                   pilot$q50_6,
                   pilot$q50_7,
                   pilot$q50_8))
## alpha = 0.86 

## Make scales! 
pilot <- pilot %>% 
  rowwise(.) %>% 
  mutate(other_trait = mean(c(other_american,
                              other_intelligent,
                              other_honest,
                              other_openminded,
                              other_generous,
                              q50_6,
                              q50_7,
                              q50_8),
                            na.rm=T),
               aff_pol = mean(c((other_therm/100),
                          (other_trait/5),
                          (q43/5)),
                        na.rm=T)) %>%  
  ungroup(.) 

## alpha for joint scale? 
psych::alpha(cbind(pilot$other_therm/100,
                   pilot$other_trait/5,
                   pilot$q43/5))
## alpha =  0.82 

## Run some models 
m.ap <- lm(aff_pol ~ as.factor(treatment), 
           data = pilot) 

## Item-by-item models 
m.ft <- lm(other_therm ~ as.factor(treatment), 
           data = pilot) 
m.trait <- lm(other_trait ~ as.factor(treatment), 
              data = pilot) 
m.trust <- lm(q43 ~ as.factor(treatment), 
              data = pilot) 

## Is one treatment more effective? 
t1 <- lht(m.ap, "as.factor(treatment)1 = as.factor(treatment)2") 
t2 <- lht(m.ft, "as.factor(treatment)1 = as.factor(treatment)2") 
t3 <- lht(m.trait, "as.factor(treatment)1 = as.factor(treatment)2") 
t4 <- lht(m.trust, "as.factor(treatment)1 = as.factor(treatment)2") 
## no difference 

## Output to a table  
models <- list()
models[["AP Index"]] <- m.ap  
models[["Out-Party FT"]] <- m.ft 
models[["Out-Party Traits"]] <- m.trait
models[["Out-Party Trust"]] <- m.trust  

## add the p-values from the test comparing the treatments & N 
rows <- tribble(~term, ~"APIndex",~"OutPartyFT",~"OutPartyTraits",~"OutPartyTrust", 
                'Iden. Treat. More Effective?','N','N','N','N',
                '(p-value)',round(t1$`Pr(>F)`[2],2),round(t2$`Pr(>F)`[2],2),round(t3$`Pr(>F)`[2],2),round(t4$`Pr(>F)`[2],2), 
                'N',length(m.ap$residuals),length(m.ft$residuals),length(m.trait$residuals),length(m.trust$residuals))
attr(rows,'position') <- c(7,8,9)

## clean up the GOF output (just N & R-Squared)
gm <- modelsummary::gof_map 
gm$omit <- TRUE 
gm$omit[5] <- FALSE 

msummary(models,
         add_rows = rows, 
         coef_map = c("(Intercept)" = "Constant",
                      "as.factor(treatment)1" = "American Identity Prime",
                      "as.factor(treatment)2" = "9/11 Prime"), 
         gof_map = gm, 
         stars = T, 
         output =  "tables/raw output/table_3_1.html") 

## Heterogeneity by Race? 
m.white <- lm(aff_pol ~ as.factor(treatment)*white, 
              data = pilot)
t1 <- lht(m.white,"as.factor(treatment)1 + as.factor(treatment)1:white = 0")
t2 <- lht(m.white,"as.factor(treatment)2 + as.factor(treatment)2:white = 0") 
## effect is for white respondents 

## Heterogeneity by Party ID? 
m.pid <- lm(aff_pol ~ as.factor(treatment)*dem, 
            data = pilot)
t3 <- lht(m.pid,"as.factor(treatment)1 + as.factor(treatment)1:dem = 0")
t4 <- lht(m.pid,"as.factor(treatment)2 + as.factor(treatment)2:dem = 0") 
## effect is for Republicans 

## Output this to a table 
## Output to a table  
models <- list()
models[["Race"]] <- m.white  
models[["Partisanship"]] <- m.pid 

## add the p-values from the test comparing the treatments & N 
rows <- tribble(~term, ~"Race",~"Partisanship",
                'Non-White/Dem Treat Effect? (Amer. Iden.)','N','N',
                '(p-value)',round(t1$`Pr(>F)`[2],2),round(t2$`Pr(>F)`[2],2),
                'Non-White/Dem Treat Effect? (9/11)','N','N',
                '(p-value)',round(t3$`Pr(>F)`[2],2),round(t4$`Pr(>F)`[2],2), 
                'N',length(m.white$residuals),length(m.pid$residuals))
attr(rows,'position') <- c(19,20,21)

## clean up the GOF output (just N & R-Squared)
gm <- modelsummary::gof_map 
gm$omit <- TRUE 
gm$omit[5] <- FALSE 

msummary(models,
         add_rows = rows, 
         coef_map = c("(Intercept)" = "Constant",
                      "as.factor(treatment)1" = "American Identity Prime",
                      "as.factor(treatment)2" = "9/11 Prime",
                      "white" = "White Respondent",
                      "as.factor(treatment)1:white" = "American Identity*White",
                      "as.factor(treatment)2:white" = "9/11 Prime*White",
                      'dem' = "Democrat",
                      "as.factor(treatment)1:dem" = "American Identity*Democrat",
                      "as.factor(treatment)2:dem" = "9/11 Prime*Democrat"), 
         gof_map = gm, 
         stars = T, 
         output =  "tables/raw output/table_three_two.html") 

###############################
## Sports Fandom Experiments ##
############################### 

sports <- read_csv(file="data/sports_experiment_data.csv")

## PID and gender (needs to be recoded)
sports <- sports %>% 
  mutate(party_id = ifelse(q37 == 1 & q19 == 1, 1,
                           ifelse(q37 == 1 & q19 == 2,2,
                                  ifelse(q37 == 2 & q21 == 1,3,
                                         ifelse(q37 == 2 & q21 == 2,5,
                                                ifelse(q37 == 3 & q20 == 2,6, 
                                                       ifelse(q37 == 3 & q20 == 1,7,NA)))))), 
         party_id = ifelse(q37 == 2 & is.na(q21),4,party_id),
         dem = ifelse(party_id<4,1,0),
         sp = ifelse(party_id == 1 | party_id == 7,1,0),
         sorted = ifelse((party_id < 4 & q38 < 4) | (party_id > 4 & q38 > 4),1,0),
         female = ifelse(q71 == 1,0,
                         ifelse(q71 == 2,1,NA))) 

## Code up the treatment conditions 
sports <- sports %>% 
  ## qt_treat is the "raw" treatment 
  ## 1 : Democrat, same team 
  ## 2 : Democrat, no team 
  ## 3 : No party, same team 
  ## 4 : Republican, same team 
  ## 5 : Republican, no team 
  mutate(qt_treat = ifelse(!is.na(q32_4),1, 
                           ifelse(!is.na(q125_4),2,
                                  ifelse(!is.na(q119_4),3,
                                         ifelse(!is.na(q122_4),4,
                                                ifelse(!is.na(q128_4),5,NA))))),
         ## this is the treatment matching the design document 
         ## 1 : same team, same party 
         ## 2 : same team, no party 
         ## 3 : same team, other party 
         ## 4 : no team, same party 
         ## 5 : no team, other party 
         treat = ifelse(qt_treat == 1 & party_id < 4, 1,
                        ifelse(qt_treat == 1 & party_id > 4, 3,
                               ifelse(qt_treat == 2 & party_id < 4, 4,
                                      ifelse(qt_treat == 2 & party_id > 4,5,
                                             ifelse(qt_treat == 3, 2, 
                                                    ifelse(qt_treat == 4 & party_id < 4,3,
                                                           ifelse(qt_treat == 4 & party_id > 4,1,
                                                                  ifelse(qt_treat == 5 & party_id < 4, 5,
                                                                         ifelse(qt_treat == 5 & party_id > 4,4,NA))))))))),
         ## dummies for treatment 
         cond1 = ifelse(treat == 1,1,0),
         cond2 = ifelse(treat == 2,1,0),
         cond3 = ifelse(treat == 3,1,0), 
         cond4 = ifelse(treat == 4,1,0),
         cond5 = ifelse(treat == 5,1,0))

## Code up the DVs 
## rs_xx is the variable rescaled to the 0-1 index 
sports <- sports %>% 
  mutate(feeling_therm = ifelse(!is.na(q32_4),q32_4, 
                                ifelse(!is.na(q125_4),q125_4,
                                       ifelse(!is.na(q119_4),q119_4,
                                              ifelse(!is.na(q122_4),q122_4,
                                                     ifelse(!is.na(q128_4),q128_4,NA))))), 
         trust = ifelse(!is.na(q24_1),(-1*q24_1)+6,
                        ifelse(!is.na(q126),(-1*q126)+6,
                               ifelse(!is.na(q120),(-1*q120)+6,
                                      ifelse(!is.na(q123),(-1*q123)+6,
                                             ifelse(!is.na(q129),(-1*q129)+6,NA))))),
         rs_ft = feeling_therm/100,
         rs_trust = trust/5,
         patriotic = ifelse(!is.na(q25_1),q25_1,
                            ifelse(!is.na(q127_1),q127_1,
                                   ifelse(!is.na(q121_1),q121_1,
                                          ifelse(!is.na(q124_1),q124_1,
                                                 ifelse(!is.na(q130_1),q130_1,NA))))), 
         intelligent = ifelse(!is.na(q25_2),q25_2,
                              ifelse(!is.na(q127_2),q127_2,
                                     ifelse(!is.na(q121_2),q121_2,
                                            ifelse(!is.na(q124_2),q124_2,
                                                   ifelse(!is.na(q130_2),q130_2,NA))))),
         honest = ifelse(!is.na(q25_3),q25_3,
                         ifelse(!is.na(q127_3),q127_3,
                                ifelse(!is.na(q121_3),q121_3,
                                       ifelse(!is.na(q124_3),q124_3,
                                              ifelse(!is.na(q130_3),q130_3,NA))))),
         open = ifelse(!is.na(q25_4),q25_4,
                       ifelse(!is.na(q127_4),q127_4,
                              ifelse(!is.na(q121_4),q121_4,
                                     ifelse(!is.na(q124_4),q124_4,
                                            ifelse(!is.na(q130_4),q130_4,NA))))),
         generous = ifelse(!is.na(q25_5),q25_5,
                           ifelse(!is.na(q127_5),q127_5,
                                  ifelse(!is.na(q121_5),q121_5,
                                         ifelse(!is.na(q124_5),q124_5,
                                                ifelse(!is.na(q130_5),q130_5,NA))))),
         hypocritical = ifelse(!is.na(q25_6),(-1*q25_6)+6,
                               ifelse(!is.na(q127_6),(-1*q127_6)+6,
                                      ifelse(!is.na(q121_6),(-1*q121_6)+6,
                                             ifelse(!is.na(q124_6),(-1*q124_6)+6,
                                                    ifelse(!is.na(q130_6),(-1*q130_6)+6,NA))))),
         selfish = ifelse(!is.na(q25_7),(-1*q25_7)+6,
                          ifelse(!is.na(q127_7),(-1*q127_7)+6,
                                 ifelse(!is.na(q121_7),(-1*q121_7)+6,
                                        ifelse(!is.na(q124_7),(-1*q124_7)+6,
                                               ifelse(!is.na(q130_7),(-1*q130_7)+6,NA))))),
         tr_mean = ifelse(!is.na(q25_8),(-1*q25_8)+6,
                          ifelse(!is.na(q127_8),(-1*q127_8)+6,
                                 ifelse(!is.na(q121_8),(-1*q121_8)+6,
                                        ifelse(!is.na(q124_8),(-1*q124_8)+6,
                                               ifelse(!is.na(q130_8),(-1*q130_8)+6,NA)))))) %>% 
  ## now find the average trait levels 
  rowwise(.) %>% 
  mutate(avg_trait = mean(c(patriotic,intelligent,honest,open,
                            generous,hypocritical,selfish,tr_mean),na.rm=T),
         rs_at = avg_trait/5, 
         affpol_index = mean(c(rs_ft,rs_trust,rs_at),na.rm=T))

## Do the trait ratings hang together? 
psych::alpha(cbind(sports$patriotic, sports$intelligent,
                   sports$honest, sports$open, 
                   sports$generous, sports$hypocritical,
                   sports$selfish,sports$tr_mean))
## alpha = 0.87 

## Do the three AP measures hang together? 
psych::alpha(cbind(sports$rs_ft,sports$rs_trust,
                   sports$rs_at))
## alpha = 0.81 

## Make the graph by condition 
pdf(file="figures/chi_levendusky_fig03008.pdf")
sports %>% 
  group_by(treat) %>%
  summarise(avg_ft = mean(feeling_therm, na.rm = T),
            sd_ft = sd(feeling_therm, na.rm = T),
            count = n(), 
            se_mean = sd_ft/sqrt(count)) %>%
  ggplot(data = .) + 
  geom_point(aes(x=avg_ft,y=treat)) + 
  geom_errorbarh(aes(y = treat, xmin = avg_ft - se_mean, xmax = avg_ft + se_mean, height = 0.001)) + 
  xlim(c(45,85)) + 
  theme_bw() + 
  ggtitle("") + 
  xlab("Average Feeling Thermometer Rating") + 
  #ylab("") + 
  scale_y_continuous(name = "", 
                     breaks = c(1:5),
                     labels = c("Same Party, \n Same Team", 
                                "No Party, \n Same Team",
                                "Other Party, \n Same Team",
                                "Same Party, \n No Team",
                                "Other Party, \n No Team")) + 
  theme(panel.grid = element_blank()) 
dev.off() 

## Run regressions for the appendix 
## First do the results for the feeling thermometers only 
m1<- lm(feeling_therm ~ cond3,
        data = sports, 
        subset = (treat == 3 | treat == 5))
m2 <- lm(feeling_therm ~ cond1 + cond3 + cond4,
         data = sports, 
         subset = (treat != 2)) 
lht(m2, c("cond1 - cond3 = cond4"))

stargazer(m1, m2, 
          column.labels = c("Model 1","Model 2"),
          covariate.labels = c("Same Party, Same Team", "Other Party, Same Team","Same Party, No Team","Constant"),
          digits = 2, 
          keep.stat = c("n","rsq"),
          type = "html",
          out = "tables/raw output/table_appendix_3_1.html")

## Now do the results for all three items 
## scale them together and form an index [0-1]
m3<- lm(affpol_index ~ cond3,
        data = sports, 
        subset = (treat == 3 | treat == 5))
m4 <- lm(affpol_index ~ cond1 + cond3 + cond4,
         data = sports, 
         subset = (treat != 2)) 
lht(m4, c("cond1 - cond3 = cond4"))

stargazer(m3, m4, 
          column.labels = c("Model 1","Model 2"),
          covariate.labels = c("Same Party, Same Team", "Other Party, Same Team","Same Party, No Team","Constant"),
          digits = 2, 
          keep.stat = c("n","rsq"),
          type = "html",
          out = "tables/raw output/table_appendix_3_2.html")
